% This function values the posterior value for the observables given the
% state estimates. 
% This is done in the case where the measurement equations reads:
% y(t) = g(x(t)) + w(t) and E(w(t)) = 0 and Var(w(t)) = Rw

function [yhatOut] = CDKF_y_posterior(yData,model,setupFilter,xhatOut,SxOut,version)

h2             = 3;                        % Squared divided-difference step size
h              = sqrt(h2);                 % Divided difference step-size
[dimy,samples] = size(yData);              % The samples length and number of observables
dimx           = size(xhatOut,1);          % # of states

% Allocating memory
yhatOut      = zeros(dimy,samples);

% Recall that the first element in x_estimated is for the steady state, x0;
for k=1:samples
    xhat         = xhatOut(:,k);
    if version == 1
        % Simply evaluate model at xhat, i.e. no account of state uncertainty
        yhat = modelFit(xhat,model);
        yhatOut(1:length(yhat),k)   = yhat;
    elseif version == 2
        % Integrate out xhat using Sx, i.e. we account for state uncertainty
        SxHat        = squeeze(SxOut(:,:,k));
        y0All        = modelFit(xhat,model);
        yhat         = ((h2-dimx)/h2)*y0All;
        for kx=1:dimx
            sypAll      = modelFit(xhat+h*SxHat(:,kx),model);
            symAll      = modelFit(xhat-h*SxHat(:,kx),model);
            yhat        = yhat +(sypAll+symAll)/(2*h2);
        end
        yhatOut(:,k)   = yhat;
    end
end

end
